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Abstract 

Although the existence of multiple stable phenotypes of living organisms enables random switching between phenotypes 
as well as non-random history dependent switching called hysteresis, only random switching has been considered in prior 
experimental and theoretical models of adaptation to variable environments. This work considers the possibility that 
hysteresis may also evolve together with random phenotype switching to maximize population growth. In addition to 
allowing the possibility that switching rates between different phenotypes may depend not only on a continuous 
environmental input variable, but also on the phenotype itself, the present work considers an opportunity cost of the 
switching events. This opportunity cost arises as a result of a lag phase experimentally observed after phenotype switching 
and stochastic behavior of the environmental input. It is shown that stochastic environmental variation results in maximal 
asymptotic growth rate when organisms display hysteresis for sufficiently slowly varying environmental input. At the same 
time, sinusoidal input does not cause evolution of memory suggesting that the connection between the lag phase, 
stochastic environmental variation and evolution of hysteresis is a result of a stochastic resonance type phenomenon. 
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Introduction 

Adaptation of organisms to time-varying and often uncertain 
environments is a classical problem in evolutionary biology. 
Existence of multiple phenotypes and random switching between 
them establishes phenotypic diversity within the population and 
has been suggested as a form of bet-hedging strategy that increases 
the chances of survival and growth rates of the total population 
[1]. It is intuitively clear that selection should favor those 
phenotype switching probabilities that match in some way 
environmental variation rates. To see what might be a good 
choice of the switching probabilities, Kussel & Leibler [2] 
compared random switching between phenotypes to responsive 
switching using a continuous time model for a discrete-valued 
environment by imposing a cost on the non-random responsive 
strategy and assuming that random switching rates are indepen- 
dent of the environment. Others [3] considered a different 
approach where switching probabilities were viewed as a single 
valued function of a binary environmental variable that may be 
favorable or unfavorable to a particular phenotype in terms of the 
phenotype growth rates. This work concluded that, under some 
circumstances, small switching probability from favorable to 
unfavorable phenotype may be advantageous for the growth of 
the entire population. In their experimental work [4-6] the same 
group was able to tune the phenotype switching probabilities 
utilizing bi-stability in the galactose utilization network of 
Saccharomyces cerevisiae obtaining agreement with the model. 



Phenotypic multi-stability in biological systems is related not 
only with bet-hedging behavior, but also with persistent memory 
of history called hysteresis. The term "hysteresis" seems to have 
been coined by James Alfred Ewing [7] in connection with the 
ability of some magnetic materials to retain their magnetization 
state long after the magnetizing magnetic field has been removed. 
Today it is used much more broadly to refer to any memory based 
relationship between an input and state of a system that does not 
depend on the rate at which the input varies in time [8] . The most 
basic, yet non-trivial, hysteresis is exemplified by a bi-stable relay 
illustrated in Figure 1, where the current state of the relay is 
determined not only by the external input, but also by the previous 
relay state. The key characteristic of the bi-stable relay hysteresis is 
the threshold separation, which is the difference between the 
values of the external input at which the state switches. 

Classical example of a bi-stable relay hysteresis in biological 
systems is the history dependent behavior of the lac-operon studied 
in E. coli bacteria. Lac-operon can be viewed as a collection of 
genes associated with transport and metabolism of lactose. Novick 
& Weiner [9] as well as of Cohn & Horibata [10-12], relying on 
prior work of others [13-17], demonstrated that two phenotypes 
each associated with "on" and "off' state of the lac-operon 
expression can be obtained from the same culture of genetically 
identical bacteria. The fraction of each corresponding sub- 
population depended on the history of the exposure to the 
inducer. Similarly to the relay illustrated in Figure 1, the lac- 
operon state was induced (switched on) when the extracellular 
inducer concentration (input) exceeded an upper threshold. The 
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Figure 1. Illustration of bi-stable relay hysteresis. When the state is 1, the relay switches to state 0 at the lower input threshold, while starting 
from state 0 the relay switches to state 1 at an upper threshold. Thus, as long as the current input is within the bi-stability range (i.e. between the 
thresholds), the relay remembers whether the input has entered this range from below or from above. Panel A shows the input-state diagram. Panel 
B presents an example of input graph. The vertical interval between the horizontal dashed lines on panel B corresponds to the horizontal bi-stability 
interval on panel A. The intersection of the lower (upper) dashed line with the input graph on B defines a moment of switching from state 1 to 0 
(from 0 to 1). 
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operon was switched off when the inducer concentration fell below 
the lower threshold. Both phenotypes remain stable through 
multiple generations of the bacterial culture after the extracellular 
concentration of the inducer is reduced to lower levels, but not 
removed completely. Novick & Weiner did not use the term 
"hysteresis" to describe their observations, but effectively that is 
what it was. 

It has been pointed out a number of times that random 
switching and hysteretic switching (memory) can both be observed 
in biological systems having multiple stable phenotypes [18-32]. 
Consider, for example, bi-stable relay of the galactose utilization 
network of Saccharomyces cerevisiae studied in [6,33] and used in 
experiments related to advantages of population diversity [34]. For 
a certain set of control parameters, the bi-stability displayed 
pronounced random switching between phenotypes. An illustra- 
tion that does not rely on concepts related to thermodynamic 
equilibrium is shown in Figure 2. The figure shows an S-shaped 
curve whose features may be randomly changing due to random 
variations in gene expression, for example. Stable phenotype states 
can be found on the horizontal segments of this curve, while the 
segment with the negative slope is unstable. If one considers a 
binary environmental variable, the random changes in the curve 
features may cause random transitions from one stable state to 
another with different probabilities at different values of the binary 
environment. When the value of the environmental variable is 
E = E\, the transition from the state on the lower part of the curve 
to the upper part occurs with greater probability. On the other 
hand, when E = Eo, transitions from the upper to the lower 
branch are more likely. Such transitions may relatively quickly 
randomize the phenotype erasing any memory of the past 
phenotype state [35-41]. 

In the above description, one may treat the transition 
probability solely as a function of the environment because the 
environment is binary valued. This becomes more obvious for a 
different set of the regulatory network parameters when the 
galactose utilization network of Saccharomyces cerevisiae exhibits 
persistent memory over long period of time [6] . Random switching 
with rates similar to the previous case would occur in this case too, 
however, only when the input is sufficiently close to one of the two 
thresholds associated with the hysteresis. Therefore, both memory 



and random switching of phenotype can be observable if a 
continuum of the environmental input states was considered, some 
far from the thresholds and some close. Hence, if one wants to 
investigate the possibility and effects of hysteretic memory, one 
should assume the dependence of the transitions probabilities on 
the phenotype as well as on the value of the continuous 
environmental variable. This is exactly the approach taken in 
the present work. Thus, while the previous work focused on the 
random phenotype switching showing that appropriate choice of 
switching probabilities is advantageous in differently varying 
binary environments, in this work we focus on the possibility that 
hysteresis and memory can provide certain advantages in 
population growth by considering continuously varying environ- 
mental input. Rephrasing, the previous analysis dealt with the 
following optimization problem: given a certain switching rate to 
the favored phenotype for every possible state of the environment, 
find the value of the lower (typically small) transition rate to the 
unfavored phenotype that would maximize the growth of the total 
population. In this work, we are solving a different optimization 
problem by allowing a range of environmental states (between the 
thresholds) within which both switching rates to and from the 
favored phenotype are small (possibly zero) and using this range as 
an optimization parameter to maximize the growth rate. 

That is, the question we ask is: how strong of a hysteresis effect 
may one expect to evolve? More precisely, what would be the most 
advantageous threshold separation and how is it affected by the 
variability and uncertainty of the environmental input? 

This problem is studied below using a hybrid system model very 
similar to one used to describe population diversification in the 
experiments with Saccharomyces cerevisiae [3,4]. However, the 
proposed model differs in two significant ways from the previous 
model. 

One difference is that the proposed model explicidy introduces 
a non-growing phenotypes in addition to the two phenotypes 
capable of some growth in any environment being considered. 
The non-growing phenotypes are introduced instead of explicidy 
imposing a time delay as a model for lag phase in the growth of 
phenotypes after each switching event. On the one hand, using a 
non-growing phenotype is a self-consistent approach as the growth 
of each phenotype is not explicitly dependent on time. On the 
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Figure 2. Illustration of random phenotype switching. The S-shaped curve consists of two horizontal segments and a slanted segment whose 
slope may be randomly changing. Stable phenotype states can be found on the horizontal segments of this curve. For a given value of the 
environmental input E, the transition from the state on the lower part of the curve to the upper part occurs when the meeting point of the slanted 
segment with the lower state line shifts to the left of E; the transition from the upper to the lower branch occurs when the meeting point of the 
slanted segment with the upper state line shifts to the right of E. When E = E\, transitions from state 0 to state 1 occur with higher probability than 
transitions from state 1 to the state 0. When E = E 0 , transitions from state 0 to state 1 are less likely than transitions from state 1 to state 0. 
doi:1 0.1 371 /journal.pone.01 03241 .g002 



other hand, this approach simplifies the mathematical treatment 
by adding two random differential equations rather than using a 
smaller number of random differential equations with time delay. 

The key difference between the previous and the proposed 
model is the possibility that switching rates between phenotypes 
depend not only on the environment, but also on the phenotype 
itself This effectively implies that environmental input values 
(thresholds) at which different phenotypes change their switching 
rates are permitted to be different in the proposed model and is 
exactly what allows hysteretic memory to exist as a possible 
solution to the growth optimization problem. This optimization 
mechanism is only possible if the environment admits values 
between the thresholds, which is part of the reason that we model 
the variations of the environment by a continuous process. 

With these modifications the growth model becomes: 



x' = y i (E)x — k\ (E)x + Sw 
y' = y 2 (E)y-k 2 (E)y + Sz 
z' = ki(E)x-Sz 
W = k 2 (E)y — Sw 



(1) 



where "prime" denotes time derivative; k\{E) and k 2 (E) are 
environment dependent rates at which organisms switch from 
phenotypes x and y, respectively, while y\{E) and y 2 {E) are the 
corresponding phenotype growth rates; z and w represent non- 
growing phenotypes; and, 1/5 is the lag phase characteristic time. 



Different functional dependencies of the coefficients in model (3) 
can be considered. Continuous functions (piecewise linear 
sigmoidal) are used here to model dependence of growth 
coefficients on the environment: 



y + a, E<0 
■;,(£') 4-<7(£-0.5) + 0.5(7+Y, 0<£<1; y 2 (E) = 2y + a-y l (E) (2) 
Y. E>\ 



where Y i- s the minimum growth rate possible and y + a is the 
maximum possible growth rate. Effectively, a is the maximal 
growth rate advantage of the alternate phenotype. 

Although the above dependence is employed here primarily to 
illustrate key features of the model, it can be viewed as reasonable 
because 1) the growth rates of both phenotypes can be the same at 
some value of the environmental variable (which is set to £" = 0.5 
here) and 2) only a partial favoring of one phenotype over the 
other is possible when the environmental variable deviates from 
£ = 0.5 by a small amount. On the other hand, at large deviations 
of the environmental variable from its average, the difference in 
favoring one phenotype over another is bounded by some value a. 

Similarly to the previously considered models, the dependence 
of the switching rates k\(E), and k 2 (E) on the environmental 
variable will be described by step functions. However, in contrast 
to the previous model, the thresholds for these steps will not be 
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required to be the same for the two growing phenotypes. We use a 
parameter a to specify the separation of thresholds: 

( k u , E<0.5 + a ( k f , E<0.5-a /s 

k\(E)={ k 2 {E)=\ ' 3 

\k f ,E>0.5 + a \k u ,E>0.5-a W 

where k u >0 is the rate at which bacteria switch from favored to 
unfavored phenotype and kf>k u is the rate of switching from 
unfavored to favored phenotype. When a = 0, the change of the 
transition rate from one phenotype to another coincides with the 
change of the phenotype growth status from favored to unfavored 
or vice versa. This change of transition rates can be characterized 
as "realistic" strategy. Positive a implies that there is an interval of 
the environmental input 0.5 — cc<E<0.5 + a over which both 
growing phenotypes have low transition rates k u . In this case one 
phenotype retains its low transition rate even after its growth status 
has changed to unfavored, while the other reduces its transition 
rate even before its growth status changes to favored. For this 
reason the strategy corresponding to positive a can be character- 
ized as optimistic. Negative a means that both phenotypes have 
high transition rate kf over the interval 0.5 + a<is<0.5 — a 
which corresponds to a pessimistic strategy, although somewhat 
different from bet-hedging. 

Variation of the environment will be modeled here in two 
distinct ways: by a random process and by a periodic symmetric 
{E(t) = 1 —E(t — T/2)) function. One of the random processes 
employed here is the Ornstein-Uhlenbeck (OU) process that 
describes diffusion-like motion in a one-dimensional parabolic 
potential centered at the point £" = 0.5 where the growth rates of 
the phenotypes are equal: 



dE=-a(E-0.5)dr + dW, 



(4) 



where a is a stiffness parameter associated with the parabolic 
potential well and dW t is the derivative of the Wiener process 
(white noise) creating stochastic fluctuations around the point 
£ = 0.5. Average time xe of passage of the interval 
0.5 — |a| <£'<0.5 + |a| by the OU process can be viewed as a 
certain characteristic time of this process. A modification of the 
OU process that corresponds to a double-well, rather than a 
parabolic potential will also be considered: 



Y=yN, Z = zN, W = wN, system (1) becomes 

X' = y 1 (E)(l - X)X - ki (E)X + 3W- y 2 (E)XY 
r = y 2 (E)(l -Y)Y-k 2 (E)Y + 3Z- y x (E)XY 
Z' = ki (E)X -SZ- y x (E)XZ - y 2 (E) YZ 
W = k 2 (E) Y- 3 W- Y] (E)XW - y 2 (E) YW 



(7) 



While populations in model (1) grow exponentially, solutions of 
model (7) fluctuate near a positive equilibrium. 

Simulations of system (7) were performed using the Runge- 
Kutta method with the input E = E(t) obtained by the Euler 
method. A typical time step was h= 10 ~ 3 ; the stochastic term in 
(4) was modeled by a sequence of independent random variables 
rj t ~ N(0,h). The time interval T = 2000 of individual simulations 
was chosen sufficiently large to ensure the convergence of the 
growth rate 



m= 7 [' (ym^m^+ymx))Y(x))dx 

t Jo 



(8) 



dE = (a(E- 0.5) - (E - 0.5) 3 )dt + dW, 



(5) 



to its asymptotic value (6). Figure 3 presents a typical plot of k{t) 
which has two parts corresponding to the first half and the second 
half of the time interval 0 < t < T. The first part includes the 
process of relaxation of X(f) to its stationary value. The second part 
is almost stationary with deviations from the stationary value 
within 1%. Since the variance of X(T) is this small, all the plots 
presented in Figures 4-6 look almost smooth. 

Each point of the plots in Figures 4-6 was obtained by averaging 
the value of X(T) over 20 simulations of system (7) with stochastic 
input (4). The stepping of a was Aa = 0.01. A similar procedure 
was used to obtain the plots in Figure 7A with the exception that 
input (4) was replaced by input (5). It has been checked that the 
plots shown in Figure 4 remain essentially the same when a 
smaller time step /i=10~ 5 ora longer time interval T = 20000 is 
used. 

Results 

The results are presented for zero transition rate k u = 0 from 
favored to unfavored phenotype. 

When the lag time is close to zero (large 3), model (1) can be 
approximated by the system 



Deterministic periodic environmental variation will be taken as 
sinusoidal having a half-period Xe- 

The effect of the threshold a on the population growth will be 
investigated along with the effects of parameters k u f and 3 using 
the Lyapunov's exponent to represent the asymptotic growth: 



X = lim 



' Tl (£(i))x(T)+y 2 (£(%(i) 
o x(x)+y(x) + w(x) + z(x) 



dx 



(6) 



Simulations 

In order to obtain the dependence of the average growth rate (6) 
on the threshold separation distance a and other parameters, we 
changed variables and considered the population of cells in each 
phenotype in terms of its fraction of the total population 
N = x+y + z + w. With the change of variables X = xN, 



x' = y t (E)x - k x (E)x + k 2 (E)y 
y = li (E)y + fc i (£> - k 2 (E)y 



(9) 



where the non-growing phenotypes have been removed. The 
results obtained from model (1) with 3»l and model (9) are 
similar. Figure 4A presents the dependence of the average growth 
rate X on a for the system with zero lag phase, which is driven by 
OU environmental process (4). The horizontal graph corresponds 
to kf = Q, the case of no flow between phenotypes. Here the 
Lyapunov's exponent X is close to the arithmetic mean y + a/2 of 
the saturated highest and lowest growth rates. The other three 
curves corresponding to positive transition rate kf from unfavored 
to favored phenotype demonstrate a clear maximum, which 
increases with kf. These curves converge to the horizontal graph 
as a increases, the reason being that for large positive cc the 
switching threshold values are so high that the environmental 
input reaches them rarely, hence little switching occurs and the 
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Figure 3. Convergence of the growth rate to the Lyapunov exponent. A typical evolution X(t) of the time average of the growth rate (8) over 
the increasing time interval from 0 to 1. For 1> 1000 the process ).(t) asymptotically approaches a stationary value (6). Deviations are within 1%. The 
plot was obtained using system (7) with OU input (4). 
doi:1 0.1 371 /journal.pone.01 03241 .g003 



system behaves almost as in the no flow case kf = 0. A different 
behaviour is observed for large in absolute value negative a. In this 
limit, the environmental input is unable to pass the threshold 
values so as to turn off the switching rate, hence each phenotype is 



constandy transitioning into the other with the effect that the 
populations are permanendy mixing at constant rate kf. The 
Lyapunov's exponent X first increases and then decreases with 
increasing kf for large in absolute value negative a. 
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Figure 4. The Lyapunov's exponent X for different a values. Plots on panels A and B were obtained for systems (9) and (1), respectively, with 
OU environmental input (4). A: The curves correspond to kj =0,0.3, 2.6, 100. For a = 0 the growth rate A increases with kf, B: The curves are for 
k/=0, 0.5, 1, 1.5. For negative a the growth rate k decreases with kf. Other parameters are k u = 0, a=\, a=l, y=0.1, <5 = 1. 
doi:1 0.1 371 /journal.pone.01 03241 .g004 
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Figure 5. The effect of changes in the parameters a and <> on the Lyapunov's exponent X. Each plot demonstrates a positive optimal value 
of the threshold parameter ot, which maximizes X for model (1). A: The effect of altering the average lag time 1/5. The curves correspond to 
<5 = 2.8, 5.3, 10.3. The growth rate X increases with S. Solid lines are plotted for kf = 2.6, a relatively slow transition rate; dashed lines are plotted for 
kf = 100, a high transition rate. Other parameters are k u = 0, <x=l, a = l and y = 0.1. B: The effect of varying the stiffness of the potential well a. The 
lines correspond to o = 2, 5, 21. The growth rate X decreases with a. Other parameters are k u = 0, <r = 4, kf = 2.6, <5= 1 and y = 0.1. 
doi:1 0.1 371 /journal.pone.01 03241 .g005 




Figure 6. Results of modifying the difference of the growth rates. The parameter a measures the difference of the growth rates of the fully 
favored and fully unfavored phenotypes. The growth rate X increases with a. The curves correspond to (7 = 2,4, 8. Other parameters are fc„=0, 
/c/ =2.6, 5 = 1, a=\ and y = 0.1. 
doi:1 0.1 371 /journal.pone.01 03241 .g006 
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Figure 7. Dependence of the Lyapunov's exponent X on the parameter a for alternative environmental inputs. A: Results for model (1 ) 
with the stochastic input (5). The curves correspond to k/ =0.2, 4.7, 9.2. For negative y. the growth rate X decreases with kf. B: Results for model (1) 
with the periodic environmental input E(i) = 0.5(1+ sin/). The curves correspond to /cy =0.45, 0.7, 1.7. Again, for negative a. the growth rate X 
decreases with kf. Other parameters are the same as in Figure 4. 
doi:1 0.1 371 /journal.pone.01 03241 .g007 



Figure 4A shows that when the switching rate kf is relatively 
low, the maximal growth rate X is achieved by the bet-hedging 
(pessimistic) strategy corresponding to a negative value of a, while 
for larger kf the optimal asymptotic growth occurs for the realistic 
strategy corresponding to a = 0. This behavior is consistent with 
the results of [3] where a. was always zero, but a positive transition 
rate k u played the role of an (alternative) bet-hedging mechanism. 
Increasing k u with simultaneously setting a negative a help 
increase the growth rate in the present model in case of a relatively 
low switching rate kf (not shown in the figure). 

Figures 4B and 5A illustrate the main finding of this work. 

Figure 4B shows that optimal asymptotic growth occurs at 
positive values of a when the lag time 1/(5 becomes non-zero. A 
positive a. corresponds to the presence of a bi-stability interval 
which ensures that the cells do not switch phenotype when the 
environmental variable is placed within the range 
0.5 — n <E(t)< 0.5 + a. Interestingly, the optimal value of a is 
nearly independent of the maximal switching rates as long as the 
lag phase delay is the same. 

The plot in Figure 4B is divided into three regions, oc< — 1 
(region 1), a>2 (region 2), and — l<a<2 (region 3). The 
horizontal line corresponding to kf = 0, the case where there is no 
transitions between phenotypes, is the same as in Figure 4A. In 
region 1, the Lyapunov's exponent k rapidly decreases with kf. 
The reason is that the environment is between the thresholds most 
of the time for this region, 0.5 — a <£(/)< 0.5 + a, hence the 
majority of bacteria are nearly always in a transition state due to 
the pessimistic strategy (negative a). The higher the transition rate 
kf, the higher the fraction of the total population that is stuck in 
the groups z and w that do not contribute to the growth of the 
system, hence lower X. In region 2, the value of a is also sufficiently 
large so that the environmental input mostly remains within the bi- 
stability interval. As k u = 0, both rates k\ i(E) are nearly always 
zero due to the optimistic strategy corresponding to positive a. in 
this region, hence the majority of bacteria are in the non-transition 
states X,y and the plots of X for all kf tend to the horizontal plot 
obtained for kf = 0 as oc increases. Central Region 3 is the most 
interesting as each plot A(of) corresponding to a non-zero value of 
kf achieves a distinct global maximum at some positive value of a, 
that is for the optimistic switching strategy. 



Now, we consider how system (1) responds to variations of 
parameters. Figure 5 illustrates dependence of positive value of a. 
needed to obtain maximum asymptotic growth on the on the lag 
time and on the "stiffness" parameter of the OU process and. 
Figure 5A shows that, as the lag time 1 jb increases, the optimal 
positive value of a which grants the maximal Lyapunov exponent 
also increases. There is a direct relationship between a. and the 
average time Tg required for the OU process to pass through the 
bi-stability interval. Hence, the exit time Xe required to optimize 
the asymptotic growth tunes with the lag time: x% increases with 
the increasing lag phase delay. This trend is also in agreement with 
Figure 4A presenting the limit case of zero lag time where the 
corresponding optimal a is zero or negative. 

Examining the plots showing the dependence X(%) in Figure 5B 
for several values of the stiffness of the potential well of the 
environmental input, a, we see that as the well becomes steeper 
and the environmental input is forced to spend more time around 
the point £ = 0.5 of equal favoring of the phenotypes, the value of 
the peak in the Lyapunov's exponent X decreases. When the well 
becomes sufficiently steep, the peak is lost and the Lyapunov's 
exponent converges to the value y + 0.5(7 of the average growth 
rate of the system with no transitions between phenotypes. 

In Figure 6, we vary the parameter c, which controls the benefit 
to the growth rate that bacteria in a favored phenotype gain over 
bacteria in the unfavored phenotype. Increasing the value of a has 
an effect similar to that of shortening the lag time by increasing 3, 
cf. Figure 5A. This result can be understood if we consider a as a 
penalty for being in the wrong phenotype when the environment 
changes. When the penalty becomes too high it is no longer worth 
delaying changing phenotype and becomes better to change with 
the environment using the realistic switching strategy, that is 
setting a = 0. 

Finally, we test system (1) with environmental inputs different 
from the OU process. Figure 7A presents data obtained for input 
(5) generated by the diffusion process in a double well potential. 
When the transitions rate kf is low, the Lyapunov's exponent is 
maximized by a = 0. However, as the transition rate kf increases, 
the peak in the Lyapunov's exponent profile /1(a) shifts to the 
region of positive a. That is, as in case of environmental input (4) 
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(see Figure 4B), the optimistic strategy grants more fitness to the 
population than the realistic strategy for non-zero lag times. 

Plots in Figure 7B were obtained for a periodic input, which 
represents a fully predictable deterministic pattern of environment 
variations. Here the graph X(a) follows a complex profile as a is 
varied from the region a < — 0.5, where the transition rates k\i(E) 
are always equal to kf, to the region a>0.5, where there are no 
transitions between the phenotypes (k\^(E) = 0). The average 
growth rate X has a local peak in the region cx<0 and the peak 
value increases with kf. For small transition rates kf, the value of X 
at this local peak is still less than the growth rate y + 0.5(7, which is 
achieved for a>0.5 by the regime without transitions. For larger 
kf, this peak becomes the global maximum, that is the maximum 
average growth rate is achieved by the bet-hedging (pessimistic) 
strategy corresponding to an <x<0. As kf increases further, the 
peak shifts towards the point n = 0. This behavior agrees with the 
results of [3]. However, we see that for the present model kf 
should be large enough to favor the bet-hedging strategy; 
otherwise, the negative effect of the switching cost dominates 
and the strategy forbidding transitions between the phenotypes 
becomes optimal. 

Importandy, in case of the periodic environmental input, 
positive values of a. do not help growth for any sufficiently high 
switching rate kf. This contrasts with our results for the stochastic 
environments. In the following discussion, we associate this 
difference with the fact that the time required for the periodic 
input to pass through the bi-stability interval does not depend on 
a. 

Discussion 

For the case of zero lag phase, the results obtained in this work 
are entirely consistent with the results obtained in [3,4] where only 
binary environmental signal was considered. The fact that 
environmental signal is not binary in this work does not appear 
to have any significant impact on the resulting conclusions. For 
high switching rate capacity kf » 1 ji£ we still conclude that no 
switching into unfavored phenotype is needed to hedge the bets 
because the population is capable of adjusting quickly to the 
environment. For sufficiendy low switching rate capacity kf the 
population is no longer capable of adjusting sufficiently quickly 
and bet hedging develops through non-zero switching probability 
into the unfavored state (the optimal switching rate k u becomes 
positive). Interestingly, negative threshold separation a, which 
effectively slows down switching into the favored phenotype in a 
certain range of the environmental input, can further increase the 
asymptotic growth rate. This could be viewed as an additional 
mechanism to tune the characteristic time of phenotype variation 
to the characteristic time of environmental variations. However, 
the authors were not able to find any experimental work 
suggesting such behavior actually occurs perhaps due to the fact 
that non-zero lag phase is nearly universal among micro- 
organisms. 

The main finding of this paper is that, as Figure 4B illustrates, 
when the lag phase delays the growth of any phenotype, the 
optimal phenotype switching strategy involves evolution of a 
positive threshold separation a. As discussed above, positive 
threshold separation a leads to hysteresis for relatively large values 
of phenotype switching rate kf or to nearly hysteretic behavior for 
smaller values of kf . 

Experimentally this hysteresis is observed when the environ- 
mental variable changes in time slowly and hence phenotype 
switching occurs at two distincdy different threshold values of the 



environmental state. On the other hand, at the environmental 
variation rate for which a particular choice of a maximizes the 
asymptotic growth rate, one would not observe hysteresis because 
short-term switching dynamics is not negligible. In fact, optimal 
choice of a corresponds to a form of stochastic resonance where 
the delay time associated with the lag phase is a fraction of the 
average period of the random phenotype switching caused by the 
random environmental input oscillation. Indeed, model (1) 
demonstrates the growth of the optimal threshold difference with 
increasing lag phase time 1 /S, which implies tuning of the average 
time interval between the switches Xe with the lag phase delay 1/8. 
This tuning suggests that measurements of hysteresis and lag phase 
can help characterize fluctuations of the natural environment. The 
growth maximizing relationship between the first passage time Xe 
and the lag phase delay time also depends on the effective 
difference in the growth rates of the two phenotypes. 

It is worth noting that deterministic environmental input does 
not lead to the same type of phenomenon, as illustrated by 
Figure 1 B where the time interval between changes of phenotype 
switching rates is independent of the threshold difference and is 
always equal to half the period. This is interesting because it 
suggests that threshold difference is only useful in the presence of 
environmental uncertainty helping the system to minimize the risk 
of changing its phenotype switching rate too often. 

The presence of hysteresis in model (1) becomes apparent when 
the typical rate 1 ji£ of the input variations is low compared to the 
inverse lag phase time S and the rate kf of switching to the favored 
phenotype but high (due to the low level of noise in the switching 
thresholds, see Figure 2) compared to the rate k u of switches to the 
unfavored phenotype. These conditions ensure that all bacteria 
switch quickly and almost simultaneously to the favored phenotype 
whenever the environmental input leaves the bi-stability interval 
0.5 — ol<E<0.5 + ol from either end, while there are almost no 
transitions to the other phenotypes when the environmental input 
is inside this interval. This is exactly the behavior described by the 
bi-stable relay illustrated in Figure 1 , which represents the simplest 
form of hysteresis and hysteretic memory. That is, for slow 
variations of the environmental input, model (1) demonstrates the 
same memory on the level of the whole population as we assumed 
in individual bacteria, the reason being that no interaction 
between organisms has been explicitly included in the model 
and the lag phase delay as well as separation of the switching 
thresholds are properties of an individual organism. 

The mechanism that promotes the separation of thresholds and 
its correlation with the rate of environmental variations and the lag 
phase time can be explained by a trade off between too much 
responsiveness to environmental variations (for small a), with the 
associated cost of often transitions, and too much inertia (for large 
large a), which leaves too many bacteria in unfavored states. The 
positive threshold difference a decreases the rate of transitions 
from less to more favored phenotype when favoring is not strong. 
Such suppression of back and forth switching agrees with some 
experimental findings [42,43] . A slight decrease in the growth rate 
due to small fluctuations of the environment from the point where 
both phenotypes are equally favored can be less dramatic than a 
drop in the growth rate due to passing through the lag phase 
induced by a switching event. 
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